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Abstract. We employ density functional theory to study in detail the crystallization 

of super-paramagnetic particles in two dimensions under the influence of an external 

magnetic field that lies perpendicular to the confining plane. The field induces non- 

^4-H fluctuating magnetic dipoles on the particles, resulting into an interparticle interaction 

that scales as the inverse cube of the distance separating them. In line with previous 
findings for long-range interactions in three spatial dimensions, we find that explicit 
2 inclusion of liquid-state structural information on the triplet correlations is crucial 

^H to yield theoretical predictions that agree quantitatively with experiment. A non- 

1-^ perturbative treatment is superior to the oft-employed functional Taylor expansions, 

^j truncated at second or third order. We go beyond the usual Gaussian parametrization 

O of the density site-orbitals by performing free minimizations with respect to both the 

shape and the normalization of the profiles, allowing for finite defect concentrations. 
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1. Introduction 

Classical density functional theory (DFT) is the method of choice to the study 
of inhomogeneous fluids [1]. Perhaps the most extreme inhomogeneities arise in a 
crystalline solid, where the density field p{r) is both periodic and shows extreme 
differences between its local values on the lattice sites and in the interstitial regions. 
DFT has been successfully applied to the problem of crystallization of a number 
of different systems [21 El SI El E], mostly in three spatial dimensions. Here, the 
most popular system is the prototype of hard spheres, for which a geometry-based 
theory [3 El IS] has proven quite successful. For soft interactions (THl HH |T2] , however, 
where one cannot assign geometrical measures to the interacting point-particles, one 
has to resort to other functionals. In particular, it has been shown [13] that 
for long-range interactions, structural information of the liquid on the pair-level is 
insufficient and triplet fluid correlations should be allowed to explicitly flow into the 
construction of the functional. Even less is known for crystallization in two spatial 
dimensions [HI [151 US E] • Here, we consider a combination of the two above-mentioned 
cases in considering long range interactions in two spatial dimensions and we study in 
detail the role played by accurate liquid-state information on triplet correlations in 
determining phase boundaries between a fluid and the coexisting crystal. 

In this paper, we study freezing of a classical two-dimensional model fluid, namely of 
a fluid of aligned dipoles directed perpendicular to the 2D-plane and repelling each other 
with a soft 1/r^ inverse-power pair potential, with the help of density functional theory 
(DFT). In Ref. [ITj we studied freezing of the the dipolar system with the modified 
weighted density approximation and its extension to third order correlation functions. 
Within this paper we will extend our previous study in several ways: We allow for a 
finite defect concentration and relax the constraint of Gaussian density peaks in the 
crystalline phase, as, e.g., suggested for hard sphere crystals in Ref. [H]. Furthermore, 
we systematically study the influence of perturbative and non-perturbative inclusion of 
higher order correlation functions of the liquid in the density functional approximation 
on the freezing transition. We employ two different approximations to the three-particle 
correlation functions, which lead to substantially different results, therefore signalling 
the importance of an accurate approximation of the latter. 

We use different approximations to the DFT — based on the famous and powerful 
approach by Ramakrishnan and Yussouff [IS], but extending on the latter in taking 
higher-order terms into account, as will be described below. The quantity to be 
approximated in the DFT of freezing is the excess Helmholtz free energy functional 
-^cx[p(^)], a unique functional of the inhomogeneous one-particle density p(r) of the 
solid [i\. The uniqueness property implies that the excess free energy can be formally 
expanded about the excess free energy of a homogeneous fluid at a uniform density p in 
terms of density difference Ap(r) = p{r) — p: 

/?Fe.[p(r)] = /3Fex(p) - V ^ / dri . . . drj^^ (r, . . . , r„; p) Ap{r,) . . . Ap(r„), (1) 

n=i ^- Jy 
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where /3 = l/^ksT) and V is the volume occupied of the system. -Fex(p) is the Helmholtz 
excess free energy and the Cq are the n-particle direct correlation functions of the fluid, 
which are well known up to second order for dipolar fluids [IT]. 

Within the theory of Ramakrishnan and Yussouff this series expansion is truncated 
at second order. We therefore refer to the theory as "second order theory" (SOT). 
Part of the reason for this truncation lies in poor knowledge about higher than second 
order correlation functions; the truncation is not well justified in the problem of freezing, 
since here Ap is not a small parameter. In particular, it has been extensively shown that 
the SOT fails to accurately predict freezing for systems interacting via long-range pair 
potentials for three-dimensional systems [131 ED] • We will show in this work, that also 
for the two-dimensional dipolar system the SOT highly underestimates the stability 
of the crystal. Therefore, several approaches have been employed to include higher 
than second-order terms in the expansion — in a perturbative ^U\ or non-perturbative 
way [IIl[T3l|2ll[22]. 

The simplest attempt to go beyond the SOT is to explicitly include the third 
order term in the expansion in equation ([I]), which we refer to as "third order theory" 
(TOT). Employing the TOT demands an approximate form of the three-particle direct 
correlation function Cg {r, r'; p) of the fluid. We will show here, that — given an accurate 
expression for Cg (r,r';p) — including this term substantially improves the predicted 
freezing temperature of the long-range 1/r^-fluid (in line with previous findings for 
long-range interactions in 3D [23])- 

A third approach to the DFT we follow here, is the Modified Weighted-Density 
Approximation (MWDA) [21j by Denton and Ashcroft which we already presented 
for the dipolar system in brief in a previous paper [T7]. This approach includes 
first and second order correlation functions of the fluid exactly (as in the SOT) and 
higher order correlation functions in a non-perturbative, implicit fashion. We find 
that the MWDA, in two dimensions, slightly shifts the freezing transition to higher 
temperature as compared to the SOT, still highly underestimating the stability of 
the solid state. In a fourth approach we employ the so called "extended modified 
weighted-density approximation" (EMA), as suggested in Refs. [131 E2]- Different 
from the MWDA, this approximation to the density functional now includes not only 
first and second, but also third order correlation functions of the fluid exactly (as 
in the TOT). Higher than third-order correlation functions are contained in a non- 
perturbative, implicit fashion, following a similar scheme as in the MWDA. For the 
dipolar system we find that this approach leads to a very accurate value of the freezing 
transition temperature, lying slightly above the one obtained from the simpler TOT. 
The two-particle correlation functions of the liquid are obtained from liquid state 
integral equation theory and from simulation. The three-particle correlation functions 
are obtained applying two approximations, both based on the two-particle correlation 
functions: The first approximation used is by Denton and Ashcroft (DA) [23], and the 
second is by Barrat, Hansen, and Pastore (BHP) [23] . 

We find that the inclusion of higher order correlation functions in a perturbative 
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(TOT) or non-perturbative (EMA) way subsequently increase the freezing transition 
temperature, thus broadening the range of the thermodynamical stabihty of the crystal. 
In fact, we find the freezing transition temperature to be in good agreement with 
experiment [26] and simulation [271 1211 129]. The importance of inclusion of third order 
correlation functions is addressed to the long-range nature of the dipole-dipole pair 
interaction. 

The rest of this work is organized as follows: In section 2 we give a brief description 
of the MWDA and of the EMA. In section 3 we apply the different approximations to 
the DFT to freezing of monodispers two-dimensional liquids. The theory is adapted to 
the dipolar system under study in section 4. In section 5 we present the resulting phase 
diagrams and different structural properties of the crystalline system, and we conclude 
in section 6. 

2. Modified weighted-density approximation and its extension to 
third-order correlation functions 

It is well known that the intrinsic Helmholtz free energy of an inhomogeneous system 
can be divided into an "ideal" and an "excess" part, 

F [p(r)] = Fid [p(r)] + F^. [pir)] . (2) 

The "ideal" term 

Fid [p(r)] = r' I drp(r) {in [p{r)A'] - l} , (3) 

is known exactly. In equation ([3]) A is the thermal de Broglie wavelength. The excess 
part can only be calculated approximately. In contrast to the SOT and TOT, within the 
MWDA and EMA the excess free energy of the inhomogeneous system is approximated 
by setting it equal to the excess free energy of a uniform liquid evaluated at a weighted 
density p, 

Fe. [p{r)] ^ FWe ip^r)] = N^p^/^) , (4) 

where superscripts denote the approximations to the DFT, MWDA (M) and EMA (E), 
respectively. A^ is the number of particles in the system and /o(p) is the excess free 
energy per particle of the liquid at the weighted density p. The latter is expressed as 

M/E 



P 



[p{r)] = ]^ / d^ dr'p{r)p{r')w {r - r'; p) 

+ — / drdr'dr"p(r)p(r')p(r'Xr-r',r-r";p) , (5) 

where the second term only appears in the EMA and not in the MWDA. The weight 
functions w{r; p) and v{r,r';p) are determined in such a way that the approximate 
functional Fex [p(^)] is exact up to second (MWDA) or third (EMA) order in density 
difference Ap(r) = {p{r) — p), i.e., up to that order equation Q and equation ([I| do 
agree. Note that the weighted density p is determined self-consistently, as it appears as 
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an argument of both weight functions. In order to obtain equahty with equation ([I| up 
to second or third order in Ap we demand the weight functions to be normahzed, i.e., 

/ drw{r;p) + — dr dr'v{r,r'; p) = 1 . (6) 

Jv ^ Jv Jv 

and to fulfill the requirements 



S'K 




-lJ2) 



= -p-^c','>(r-r';p), 

= -p-^cf{r-r\r-r"-p), (7) 

where Cg (r; p) and Cq (r, r'; p) are the two- and three-particle correlation functions of 
the liquid with density p which are an input to the theory. These conditions uniquely 
determine the weight functions. In order to obtain the simple algebraic equations for 
V and w that can be found in Ref. [22] a further approximation has to be made: The 
inner integral in the second term of equation (|6| is assumed to be equal to a constant, C 
(demanding the first term in equation ^ to be equal to 1 — C), where C is independent 
of the fixed space coordinate of the weight function v{r,r'; p). The weighted density p 
in equation ([s]) is independent of the choice of C [22] . 

For non-zero wave vectors {k ^ 0,k' ^ 0, or k -\- k' ^ 0), the Fourier transforms 
of the weight functions w{k; p) and v{k, k'; p) are simply proportional to the Fourier 
transforms of the second- and third-order direct correlation functions Cq {k; p) and 
c}) {k,k'; p), respectively: 

-(3-¥^\k;p) =2ap)w{k-p), 

-p-'c^^\k,k';p) = 6f;,{p)v{k,k';p). (8) 

Furthermore, equation ([s]), together with equations ^ and ([T]) guarantee fulfillment of 
the sum rules 

c^^^ik, k' = 0; p) = ^r /^(fc, -k- p) = ^^^ , (9) 

where the former is the compressibility sum rule, and where the superscripts on the 
correlation functions indicate that these functions are the Fourier transforms of the 
functional derivatives of the approximate excess free energy functionals in the limit of 
constant average density p [c.f. equation ([T])]. The primes on the excess free energy 
density /o denote derivatives with respect to density. 

Due to the self-consistency requirement, the approximate excess free energies of 
both the MWDA and the EMA include contributions from arbitrarily many higher 
orders. However, if expanded about the excess free energy of a fluid with the same 
average density as the inhomogeneous system according to equation (jT|, the MWDA 
only gives even order terms and estimates the odd order terms zero. Contrary, the 
EMA includes, approximately, contributions from all higher order terms. In particular, 
it includes the exact third-order term, which is an input to the theory. 
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3. Application of the different approximations to the DFT to freezing of 
monodisperse two-dimensional liquids 

In order to find the equilibrium one-particle density Peq(^) of a system at a given average 
density p and temperature T we minimize the approximate total free energy functional 
F[p(r)] of equation ^ with respect to the inhomogeneous one-particle density p{r) 
for fixed p. As described, for example, in Refs. [211 122] this minimization is pursued 
in a number of subsequent steps, depending on the kind of approximation: For all 
approximations to the DFT, first, an appropriate parametrization for the inhomogeneous 



one-particle density is made (we will employ a free minimization in section 5.3). Within 
the SOT and TOT, we can now, in a second step, calculate the excess and ideal parts 
of the Helmholtz free energy according to equations ([I]) and (tsl). However, within the 
MWDA and EMA, the excess part is given by equation Q, with the weighted density p 
obtained in an intermediate step according to equation ([s]). In a final step, minimization 
is carried out with respect to all free variables in the parametrization of p{r). 

The crystalline one-particle density which we expect to be in equilibrium for 
low temperature and/or high density has the symmetry of the triangular crystal — 
the quadratic lattice is thermodynamically unstable for the whole range of accessible 
densities/coupling constants and we expect mechanical instability with respect to the 
triangular lattice for any coupling. We can therefore express p{r) as a sum over 
reciprocal lattice vectors (RLV's) of the triangular lattice: 



p{r) = p 



1 + J2 l^KC^""' 



(10) 

where p is the average density of the solid, {K} is the set of reciprocal lattice vectors 
(RLV's), and where the pk are the dimensionless Fourier components. In terms of 
Fourier components the excess part to the Helmholtz free energy within SOT and TOT 
now reads 

PF!np{r)]/N = (3h{p) - ^ 5^ pVcf{k- p) 

2 

~"6"5Z 5Z f^KPK'P-{K+K')c'i\K,K';p), (11) 

the superscript referring to the SOT (S) and to the TOT (T), respectively. The third 
term only appears in the TOT. 

Within the MWDA and EMA, the weighted density, equation d5|, now reads 

] 1 + 5Z ^^^ (^; P) + P 5^ Yl f^KPK'P~{K+K') 'j^ — > 

As in equation rtsl) the three-particle term only appears in the EMA. 

Since a free minimization of the approximate Helmholtz free energy with respect 
to an infinite number of Fourier components pk at all RLV's is intractable, we make 



^M/E ^ ^ 



■ (12) 
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a simple ansatz for the one-particle density which is a superposition of normalized 
Gaussians centered around the lattice sites of the triangular lattice: 

p(r) = ^ V] exp [-a \r - Rf] , (13) 

R 

where a is the localization strength, n^ is the average number of particles occupying a 
lattice site, yielding a vacancy concentration n„ = 1 — Uc, and {R} is the set of Bravais 
lattice vectors of the triangular lattice with lattice constant a = (-\/3nc/2p)^/^. Thus, 
the Fourier components fiK now simply read 

/ix = e-""'/'- . (14) 



The ansatz, equation (13), was chosen in such a way that the system forms a triangular 
lattice for any finite a keeping its average density p fixed. For a — ;> the density profile 
becomes flat and the system turns into a liquid. We thus end up with two minimization 
parameters a and n^. 

This ansatz disregards a possible partition of the system into coexisting liquid and 
crystal phases of different densities keeping the overall average density fixed. However, 
this is accounted for by performing a common-tangent construction to the crystal and 



liquid volume free energy densities in the end. Furthermore, equation ( 13 ) disregards the 
spatial anisotropy of the density site profile at each lattice site. We wfll see in section 5.3 
where we relax the constraint on the density peaks, that both, the assumption of isotropy 
and the Gaussian shape are well justified close to the positions of the Bravais lattice 
vectors, i.e., where the density is reasonably large {p{r) > p). 



Employing the ansatz of equation (13) for the inhomogeneous density, the ideal 



part of the Helmholtz free energy [equation (^] can now be written as a function of a 
and Uc only: Fid [p(^)] = F'i4{a,nc] p). For ric = 1 it reads 

— Fid(a, Uc = 1; p) = const + In(pL^) + (?(«*), (15) 



G{a*)= / da. ^("'"^'"- = ^^ n 
'Ai P 



p{x,a*,nc = 1) 



(16) 



P 

where const is an irrelevant constant and L is a density-independent length scale of 
the system, x = rp^^'^ and a* = a/ p are the dimensionless space coordinate and 
localization strength, respectively, and the integral is performed over the area Ai of a 
unit cell. The function G{a*) is approximated for small and large localization strengths 
by its analytically known asymptotics 

^^~\G2ia*)=\n{a*/7i)-l, «* > 1 , ^^ 

where K^: = K / p are the dimensionless RLV's. For intermediate values of 2 < a* < 50 
the function G{a*) was calculated numerically. The function G{a*) and the asymptotics 
of equation ( [l7| ) are plotted as a function of a* in figure [1} 

The ideal free energy for values rzc 7^ 1 is obtained via the simple scaling relation 

— Fid(a, ric, p) = const + In(pL^) + G{nca*) . (18) 
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Figure 1. The function G{a/p) and its analytically known asymptotics for small and 
large localization strength. 



4. The dipolar system 

We now turn to the system of monodisperse particles which repel each other with an 
inverse-power pair potential u{r) = Uq/t^, where Uq is a parameter with dimensions of 
energy x volume. For the specific realization of two-dimensional paramagnetic colloids 
of susceptibility x exposed to a magnetic field B which is directed perpendicular to 
the 2D plane, we have uq = (xS)^/2 in Gaussian units [30]. Here, we assume perfect 
alignment of the magnetic dipoles with the external field which is well justified for 
xB"^ ^ ksT [31]. The thermodynamics and structure depend, due to simple scaling, 
only on one relevant dimensionless coupling parameter [32] 

3/2 

r = ^^^. (19) 



UqP 
knT 



Therefore, it is convenient to express all quantities in terms of F and consider coupling 
parameters rather than densities via this scaling relation. Correspondingly, the excess 
free energy within the SOT and TOT [equation [I] now read 

/?i^er (r, a) IN = /?/o(F) - \ Y. e-^^'^'-c^oKK- F) 



K=/=0 



6 



^-iK^+K-MK+Kr)/4a^i^)^K^ K'. T) 



pco 



and 



(20) 

Ky^O K'^O-K 

the third term only appearing in the TOT. Here, F is the coupling constant 
corresponding to the average density p according to equation (19), r — '^^ 
Cq = p^Cq are the dimensionless correlation functions of the fluid in ±v;>^ip±vj 
respectively. 

For the MWDA and EMA, the weighted coupling constants F now read 

1 



F(F,a) = F 



1 



— ^Ve 
3/3F/^(F) ^, 



-K^/2ap,{2) 



{K;V) 
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9pr /o(r) _r:-^o Kvo-x 



(K-- 



'+^"+^^')/2"c(')(K,K';f) 



3/2 



(21) 



where /o(r) is the derivative of the excess free energy density with respect to couphng 



constant. As in equation (12) the third term only appears in the EMA. 



In order to solve equations (20) and (21) we need the two- and three-particle 
correlation functions Cq {k;T) and Cq (fc; F) and the excess free energy density /o(r) 
of the corresponding liquid for a wide range of coupling constants F. The two-particle 
correlation function is obtained with liquid state integral equation theory or from 
computer simulations. In the first case, following the procedure described in Ref. 
we solve the Ornstein-Zernicke (OZ) equation 



h{k) 



s(2) 



(k) 



s(2) 



(k) 



(22) 



which relates the dimensionless Fourier transform h{k) = ph{k) of the total correlation 
function h{r) to the direct pair correlation function Cg yk), numerically. Note that the 
density has been absorbed in both the Fourier transform of the total correlation function 
h{k) and in the direct correlation function Cg (k). The total correlation function is 
connected to the pair distribution function via g{r) = h{r) + 1. 

The solution of equation (22J) for the two unknown quantities h(k) and Cg (k) 
demands a constitutive equation, the so called closure relation which for any non- 
trivial case can only be determined approximatively. Two approaches which proved 
successful for the description of fluids with long-range interactions will be applied here, 
the hypernetted chain (HNC) j34j and the Rogers- Young (RY) closure relation [35] . 
They can both be written as 

h{r) = e-^"('-) {1 + /(r)-i (e^M/W _ i) | _ i ^ (23) 

where x(^) = h(r) — Cg (r) is the indirect correlation function, /(r) = 1 — e"'*'" is 
a 'mixing function' with an adjustable parameter < ^ < oo which is either sent 
to infinity (HNC) — which is equivalent to letting /(r) -^ 1 — or chosen to guarantee 
thermodynamic consistency between virial and compressibility route to the free energy 
(RY). 



The coupled equations (22) and (23) are iteratively solved by applying the method 



of fast Fourier transforms for radially symmetric two-dimensional problems as suggested 
by Caillol et al [36] and as also summarized in appendix A of [33]. In order to reach 
rapid convergence an iteration procedure for the indirect correlation function x(r) is 
used, since its Fourier transform, x{k), decays more rapidly with increasing k than h{k). 

(2) 

The iteration scheme now consists of making an ansatz for Cg , calculating x according 



(2) 



to equation (22j), obtaining the next estimate of Cq via equation (23), inserting this 



into equation (22), etc., until convergence is obtained. 



s(2)/ 



Applying this procedure we are able to calculate Cg^ (A;; F) for coupling constants 



F much larger than the experimentally known coupling of freezing F 



/ 



10 
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Figure 2. The bond-orientational order parameter g^ir) for different coupling 
constants T as obtained by computer simulation, qq decays exponentially for coupling 
constants F < 11 indicating the system to be in the fluid state. 



which enables us to calculate the Helmholtz free energy of the system deep inside the 
thermodynamically stable crystalline region. 

More accurate pair correlation functions can be obtained from computer 
simulations. We have performed extensive Monte Carlo computer simulations [HZ] in 
a quadratic simulation box of size L x L comprising 900 particles employing periodic 
boundary conditions in order to measure the pair distribution function gs{r) = hs{r) + l, 
the subscript 's' denoting the simulation result. Since the accessible range of hs{r) is 
limited to distances r smaller than a cutoff radius re < L/2 we employed an extrapolation 
technique suggested by Verlet [3H] to obtain the complete pair correlation function: 
Verlet defined a closure relation to the Ornstein-Zernicke equation 



h{r) = hs{r) 



c'i\r) 



^YWCv) 



r < Tr 



r > r,. 



(24) 



where chnc(^) is given in equation (23). The Verlet closure relation [equation (24)] 



together with the Ornstein-Zernicke equation [equation (22)] uniquely specify the direct 



.(2) 



correlation function Cq (r) for all radii r and thus also yield the correlation function in 
reciprocal space Cq {k). As for the HNC and the RY closures the Ornstein-Zernicke 
equation and the Verlet closure were solved iteratively via the indirect correlation 
function x- Furthermore, Tc was chosen the largest root of h{r) still smaller than L/2. 
For the Verlet data we checked that the simulated system does not crystallize for 
coupling constants F < 11. Here, the freezing-criterion was chosen a non-exponential 
decay of the bond-orientational order parameter gQ{r) = {exp[i6[9{r) — 9{r')]]), where 
6{r) is the angle of the bond connecting two neighboring particles according to the 
Voronoi construction (see figure [2]). The application of the Verlet closure within the 
DFT formalism was thus restricted to the range < F < 11. 

(2) 

The Fourier transforms Cq (k) of the two-particle direct correlation functions 
obtained from the three different closure relations (HNC, RY, Verlet) are shown in 
figure [3] for F = 9, which is close to the experimentally determined coupling constant 
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Figure 3. The dimensionless Fourier transform Cg (k) of the two-particle direct 
correlation function at F = 9, plotted against kj p^l"^ . Shown are simulation data using 
the Verlet closure, liquid integral equation theory using the RY closure, and liquid 
integral equation theory using the HNC closure. The arrows indicate the positions of 
the first four reciprocal lattice vectors of the triangular lattice. 



of freezing Y f ~ 10 [26j. The HNC closure underestimates the pair structure strongly 
while the RY closure is closer to the simulation data. We also show the positions of 
the first four reciprocal lattice vectors of the triangular lattice with lattice constant 
a = (v^/2p)^/^. The value of Cq at these lattice vectors crucially influences the solid 



free energies, as can be seen from equations (20) and (21). 



Within the RY-approach the excess free energy density /o is obtained by integrating 
the compressibility which is inversely proportional to the static structure factor: 



/9/o(r) 



2 r^dr 

where the pressure P is given by 

/3P _ 2 r^ dr 
~ 3 



(5P_ 
P 



P 



p/l/3 



,(2) 



\-(^^\k = ^-V') 



(25) 



(26) 



In order to obtain the excess free energy density from the simulation data we make use 
of the relation |39] 



/?(Mex) = /? 



a/3/o ^ddU 



(27) 



d(5 dV 

between the average excess energy density (wex) = \^ii=3'^i,j) and /o and integrate the 
former. Note that for both of our approaches, the RY and the Verlet closure, the virial 
and the compressibility route are equivalent. As the energy dominates the free energy 
in the strong coupling limit, T > 1, the excess free energy density scales roughly linearly 
with coupling constant, as can be seen from figure 111 

For the EMA we need the three-particle correlation function Cg (fc,fc';p) of the 
underlying fluid for a wide range of coupling constants. We use here two conceptually 
different approximations: The first approximation is by Denton and Ashcroft [23] (DA) 
which is based on a weighted density approximation to the first order direct correlation 
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Figure 4. The excess free energy density /3/o(r) as a function of coupling constant F 
using the Verlet closure, the RY closure and the HNC closure. 



function c^^\r; p(r)) of an inhomogeneous system. The DA approach leads to an 
analytic expression of Cq in terms of the one- and two-particle correlation functions 



.(1) ;;(2) 



-^0 ) 



Cg of the liquid and their derivatives with respect to density: 



~(3)DA, 



c, 







k.k') 



r^ {\k\, \k'\) + f^^ {\k\, \k + k'\) + f^^ i\k'\, \k + k'\) 



where 
f^\k, k') 



.(1)' 



;:(2)n.^s(2)' 



4^^(fc)cr(A:') + 5r(A:)4^^(^' 



~(2),t,.~{2)i 



Jl)" 



.(1)' 



-rc'i\k)c^^\k'). 



Here, primes denote derivatives with respect to density, as above, 
approximation — by construction — fulfills the symmetry condition 

4^)°^(/c, k') = g(^)°^(fc, k + k') = gf °^(fc', k + k'). 



(28) 



(29) 



The DA 



(30) 



The derivatives Cq (k) were obtained by applying a simple finite difference method 
bearing in mind that 



2~(2)/ 



ik;p) 



sr 



dcl;>{kp-^/^;T) 



— kp 



^1/2 



dc\;>ikp-^/^;T) 
dkp-y^ 



,(2) 



-2c\,'>{kp^y'^;T) 



(31) 



~(3)DA 



We calculated Cq (fc, k') taking the direct correlation function from both the RY and 



the Verlet closure. As pointed out in Refs. [22l|2llll0] the DA model, although itself not 
derived from a free energy functional but from an approximate one-particle correlation 
function, is very similar to different approaches, all based on taking three successive 
functional derivatives of approximate free energy functionals. 

(3) 

We also employed another approximation for Cq , namely a factorization ansatz of 
Barrat, Hansen and Pastore (BHP) |23]. The approximation reads 



(3) ( I 



t{r)t{r-')t{\r -r'\) 



(32) 
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The function t(r) can be uniquely determined from the second of the sum rules in 
equation ^ which in r-space now reads 

I dr'cf{ry-p)= f dr't{r)t{r')t{\r - r'\) = ^""^ i^; p) ^ ^^^^ 



We solved equation (33) numerically for t{r) applying the method of 'steepest descent' 
as outlined in appendix B of reference pSj. As opposed to the simple finite difference 

(2)/ 

approach above the derivatives Cq (k) were now obtained by iteratively solving the 
coupled differentiated Ornstein-Zernicke equation and the differentiated RY closure 
relation, as outlined in appendix B of |23] . Since it proved difficult to reach convergence 
of the iteration procedure we did not pursue this method using the Verlet closure. 
The triplet-correlation function was then obtained by a double Fourier transform of 



equation (33) using a standard expansion in Legendre polynomials, as outlined in 



appendix A of 



In the single summation in equation (21) we consider all RLV's of absolute value 



\K\ < 33 l-f^il, where Ki is the smallest RLV of the triangular lattice — this comprises 
the first 299 stars of RLV's, which is by far sufficient to reach convergence of the single 
summation. The double summation is performed over sets of equivalent triangles of 
RLV's which are each characterized by the absolute values of the two RLV's K and K', 
and by the absolute value of there included angle. For the DA model and for the BHP 
model we include 42 sets of triangles of RLV's, where that RLV of the three RLV's, 
K, K', K — K', with the largest absolute value satisfies \K\ < 4|l^i|, which also 
guarantees convergence of the double sum. 

5. Results 

We first study the infiuence of the explicit inclusion of the triplet correlation functions 
obtained with the DA model and with the BHP model on the approximate excess 
free energy according to the TOT as compared to the simpler SOT, and according 
to the EMA as compared to the MWDA, respectively. For all six approaches we 



use the two different closure relations of Rogers and Young [equation (23)], and of 



Verlet [equation (24)], respectively. 



5.1. Gaussian profiles, no vacancies 

In order to keep things simple in the beginning we keep the number of particles occupying 



a lattice site, nc, in equation (13) fixed (i.e., Uc = 1) and thus end up with a single order 
parameter, the dimensionless localization strength a* = a/ p. 

In figure [5]^a,b) we show the weighted coupling constant and the associated excess 
free energy difference per particle between the solid and the liquid state Fex{a*)/N — fo, 
according to equation (|4]), as functions of localization strength a* for a value of F = 9 
which is close to the experimentally known value of freezing, Tf ~ 10 [26], for the 
MWDA and for the EMA, using the RY or the Verlet approach to the direct correlation 
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Figure 5. (a) The weighted couphng constant as a function of a* within the MWDA 

(2) 

and within the EMA using Cq from the RY and from the Verlet closure, and using 
Cq from the DA and the BHP model for F = 9. (b) The approximate excess free 
energy difference per particle /o(r(a*)) — /o(r) as a function of a* for the same 
approximations as in (a), (c) The approximate excess free energy difference per particle 
Fex{a*)/N - /o(a* = 0) as a function of a* obtained within the SOT and TOT using 
the same approximations for the two- and three-particle correlation functions as in 
(a,b) for r = 9. (d) Comparison of F^x obtained within the four different approximate 

(2) 

theories MWDA, EMA, SOT, and TOT using Cq ' from the Verlet closure, and using 



„(3) 



from the DA model for T = 9. 



function and using the two different approaches for the triplet correlation function, the 
DA and the BHP model. The latter are both based on the direct correlation functions 
used for the respective two-particle term. In figure [sFc) the excess free energy for the 
simpler SOT and TOT are plotted as a function of a* for the same approximations 
to the correlation functions. In figure |5](d) the non-perturbative and the perturbative 
approaches are compared. Different interesting features of the different approximations 
are observed: 

(i) For all approaches used except for those where Cq is obtained within the 
BHP model the excess free energy decreases monotonically with increasing localization 
strength a*, reaching a plateau for a* ~ 400 [c.f. figure [5|b,c)]. However, employing the 
BHP model to the triplet-correlation function leads to an increase of r(a*) and Fex(a*) 
for values of a* > 80. The former behavior is intuitively expected and has also been 
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Figure 6. The total free energy difference per particle AF(r)/A'' as a function of a* 
within the MWDA and EMA (a), and within the SOT and TOT (b) using c|,^^ from 

(3) 

the RY and from the Verlet closure, and using Cq from the DA model for F = 9. 



observed in the original MWDA [21] — localization is favoured by the excess part of the 
free energy. Once the density peaks become very narrow, a further increase of a* does 
not change the excess free energy further. On the other hand, the rise of T and of Fex 
within the BHP model is regarded as unphysical. We therefore do not consider the BHP 
model any further. 

(ii) Both within the DA model and within the BHP model (for a* < 80) the sign of 
the third term in equation (21 ) is negative, i.e., the value of T is decreased as compared 



to the pure MWDA and thus freezing is favoured. It is also interesting to note, that 



within the DA model the triplet-part in equation (21 ) is much smaller than the second 



order term while it is significantly larger within the BHP model. This same behavior 
had already been found for hard spheres in three dimensions |22j . 

(iii) Although the direct correlation functions using the RY- and the Verlet- closures 
do not differ by more than ~ 10% at the position of the most important first RLV (cf. 
figure k| the difference in T between the two is quite pronounced which is due to the 
self-consistency relation in equation (21). Furthermore, as shown in figure [stb) the 
difference in excess free energy is even more enhanced. 

(iv) Inclusion of higher than second-order terms in a non-perturbative way within 
the MWDA reduces the excess free energy as compared to the simpler SOT [cf. 
figure IsVd)]. However, inclusion of higher than third-order terms within the EMA 
increases the excess free energy with respect to the TOT. 

The total Helmholtz free energy per particle is obtained by adding to the excess 
part Fex the ideal part Fjd according to equation ^. The free energy difference per 
particle AF/N = [F{a*; T) - F{a* = 0; T)]/N is plotted in figure gas a function of a*, 
for the same value of F = 9 as in figure [s] for the SOT and TOT [figure [6]^a)] , and for the 






(2) 



MWDA and EMA [figure 6[b)], respectively, using Cq from the RY and from the Verlet 
closure, and using Cq obtained within the DA model. It is found that the different 
curves of AF/N show qualitatively very different behavior for the coupling of F = 9: 
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Figure 7. The total free energy difference per particle Ai^(r)/A'' as a function of a* 



„(2) 



,(3) 



within the EM A using Cq from the Verlet closure, and using Cq from the DA model 
for r = 9, 9.4, 9.8. 



While the free energy increases monotonically with a* within the SOT and MWDA and 
within the TOT and EMA employing the RY closure it displays a local minimum with 
respect to a* at a finite value of a* within the EMA and TOT, employing the Verlet 
closure, this local minimum even turning the deep global minimum within the TOT at 
a* ~ 213. The appearance of a global minimum at a finite value of a* corresponds 
to a thermodynamically stable crystalline state while the global minimum at a* = 
indicates a stable fiuid system. 

In figure [7], we display the total free energy obtained within the EMA employing the 
DA model with the Verlet closure for three different values of F = 9.0, 9.4, 9.8. We thus 
conclude from figure W\ — this has already been presented in Ref. [17] — that the EMA 
employing the Verlet closure and the DA model yields a transition from the fiuid to the 
solid close to F = 9.4: while for F = 9.0 the fiuid is stable as indicated by the minimal 
value at a* = 0, fiuid-solid coexistence is achieved at F = 9.4 (see the two equal minima 
in figure ItI). The solid phase, on the other hand, is clearly stable for F = 9.8. The 
localization parameter at coexistence is roughly aj^j^ ~ 100. 

The curves always displays a local minimum with respect to a* at a* = 0. This is in 
accordance with the mean-field character of any approximation to the DFT, which ignore 
fiuctuations leading to a breakdown of long-range order in one and two dimensions. 
Therefore, a first-order transition between fiuid and solid state is always predicted, i.e., 
the liquid system always has to overcome a free energy barrier in order to reach the 
thermodynamically stable crystalline state. 

The freezing and melting transition constants for the first-order phase transition 
predicted by the different approximations to the DFT, F^ and Tf, respectively, are 
obtained by using Maxwell's double tangent construction to the fluid and crystal volume 
free energy densities T'^^^F/N oc F/V, where F denotes the minimum free energy with 
respect to a, and F^/^ is proportional to the average density p of the system [c.f. 
equation (19)]. F^ and Tf correspond to the freezing and melting densities, ps and 
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Table 1. Freezing and melting parameters Fy and Fs, the widths of the coexistence 
regions AF = Fs — F/, the relative displacement parameters 7, and the pressures P at 
coexistence obtained within: the SOT with the RY closure (first row); the TOT with 
the RY closure (second row); the TOT with the Verlet closure (third row); the MWDA 
with the RY closure (forth row); the EMA with the RY closure (fifth row); the EMA 
with the Verlet closure (sixth row) , where all three-particle correlation functions were 
obtained with the DA model using the respective pair-correlation function as input. 
The last row displays experimental parameters for the isotropic-hexatic transition, 
the hexatic-crystal transition and the Lindcmann parameter, obtained from real-space 
microscopy measurements of magnetic colloids confined to an air- water interface. 
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Pj, respectively. The volume free energy density is exemplarily shown for the EMA 
using the Verlet closure and the DA model in figure [8j Within this approximation we 
obtain freezing and melting with a narrow coexistence gap AF = F^ — Fj. Table I 
summarizes the freezing/melting parameters for all the approximations made. The 
data are compared against experimental results obtained from real-space microscopy 
measurements of magnetic colloids confined to an air-water interface. The experiments 
give freezing with an intermediate hexatic phase. The liquid-solid transition has also 
been studied using numerical simulation p8l 129] yielding a slightly higher inverse 
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transition temperature between 12.0 and 12.25 but these investigations suffer from finite 
size effects. 

As becomes evident from Table I, the SOT, TOT, and MWDA are not 
quantitatively satisfying theories as they either over- or underestimate the freezing 
coupling. Note that the overestimation of the freezing coupling within SOT and MWDA 
are the reason why it is not possible to feed the "exact" pair structure into these theories. 
At such high coupling, no fluid pair structures are available since the fluid spontaneously 
crystallizes in the simulation. The EMA, on the other hand, yields results in close 
agreement with experimental data. The TOT obviously underestimates the freezing 
coupling by a factor of ~ 2. 

More detailed, structural information can be extracted from the localization 
parameter of the coexisting solid. For all approximations used we find localization 
parameters at freezing in the range 99 < a^^^iTs) < 115. Strictly speaking, the 
localization parameter has no counterpart in "real" 2D systems since the particles 
are not localized due to long range fluctuations. However, if one relates the particle 
displacements to that of their nearest neighbor, one can define a finite quantity as 
7 = p {{ui — ■Uj+i) ), where Ui and -Uj+i are the displacement vectors of neighboring 
lattice sites. Disregarding nearest-neighbor correlations {ui ■ itj+i), 7 can be estimated. 
Since the nearest-neighbor correlations (wj ■ Wj+i) are expected to be positive: 

7<2p(i.^)^2/<,. (34) 

By this relation, the localization parameter of the coexisting solid gives a prediction for 7 
which is included in Table I. From experiments, 7 is known to be close to = 0.038 [2B] . 
This was shown to be in accordance with harmonic lattice theory [H]. The EMA 
yields 7 < 0.020, i.e. the EMA roughly overestimates the localization of the particles 
by a factor of 2. 7 is smaller than the experimental value, contrarily to what was 



expected from the inequality (34). This shows that there is still a need to improve the 
theories in order to correctly predict localization properties. A similar overestimation 
of the localization is also common in weighted density approximations in three spatial 
dimensions |21j . 

Another quantity of interest, which is directly connected to the Helmholtz free 
energy is the pressure at coexistence which is also included in Table I. It is obtained via 



the equations (25), (26), (27), depending on whether the RY closure or the simulation 
data were used. 

5.2. Gaussian profiles, allowing for vacancies 

In this subsection, we relax the constraint of zero vacancy concentration, 1 — rZc = 0, 



in equation (13) and instead minimize the total free energy with respect to the two 
parameters a and ric, respectively. However, instead of calculating the phase diagram 
for all approximations to the DFT and to the pair- and triplet-correlation functions, we 
focus here on the two non-perturbative approaches, the MWDA using the RY closure 
and the EMA using the Verlet closure and the DA model. In figure [9] we plot the 



Crystallization of magnetic dipolar monolayers 



19 



>; 

-0.002 ° 
-0.004 B 

-0.006 or 

-0.008 J^ 
-0.01 a, 
-0.012 t. 




■0.99 



p[F(a) - F(a=0)]/N 




Figure 9. The total free energy difference per particle AF(r)/N as a function of a* 
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and Tic within the EMA using Cq from the Verlet closure, and using Cq from the DA 
model for F = 9.49. The left panel displays a zoom-in of the right panel. 



approximate total free energy per particle of the EMA as a function of a and n^ for 
the freezing coupling constant obtained at fixed ric = 1, T = 9.49. The minimum of 
the total free energy is slightly shifted in Uc and a from (n^ ~ 1, a* ~ 98.7) towards 



(jic ~ 0.998, a* ~ 100.5). As can be seen in figure |9](a), the difference in total free 
energy per particle between the two configurations is only of the order lO'^ksT, which 
has no influence on the phase diagram within the accuracy given in Table I. 

For the simpler MWDA, however, the vacancy concentration is substantially larger, 
which has pronounced effects on the phase diagram. In particular, we find the coupling 
constants of freezing and melting reduced to (Tf ^ 37.35, Ts ~ 37.45), the liquid being 
in coexistence with the triangular crystal at the parameters Uc ~ 0.966, a* ~ 200.5, i.e., 
the relaxation of ric improves the prediction of the freezing coupling while the Lindemann 
parameter 7 ^ 0.01 is by a factor of ~ 2/3 smaller than predicted within the simpler 
theory keeping Uc = 1 fixed which — compared to the experiment — is worse than the 
result from the constrained theory. 



5.3. Free minimization 

In this final subsection we completely remove the constraint of Gaussian density peaks. 
Instead, we minimize the density functional with respect to a free, periodic density 
field p{x,y), which has the periodicity of the hexagonal lattice with lattice constant 
a = {^/3nc/2pY^'^, as above. As laid out in Ref. |12], we minimize the density functional 
of the SOT with the RY closure with respect to p{r) by calculating the overdamped 
relaxation dynamics of a highly ordered hexagonal crystal with the help of dynamical 
DFT 031 m iSl He] according to 



dp{f, t) 



j3DV ■ \p{r,t)V 



5FW,t)] 



(35) 



dt ^ y-K , , Sp{f,t) 

where (3D is the mobility coefficient, which sets the Brownian time scale tb = (pD)^^. 
Since in this work we are only interested in the equilibrium state reached after long 
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Figure 10. Height of the density peak /5(r,r) and difference in free energy Af (F) 
as a function of F obtained from dynamical DFT using the SOT with the RY closure. 



time, tb is irrelevant in the following considerations, i.e., we use equation (35) just as 



a minimization procedure to the static DFT. Starting from an initial density profile 
p{r,t = 0), equation (35) is solved numerically for times {t/rB) ^ 10 applying a 
finite difference method and keeping the coupling constant F fixed. The maximum 
time is chosen large enough to guarantee convergence towards a (local) minimum of 
the free energy landscape. The rectangular periodic box of size L^ x Ly = ^/3a x a 
with a discretization of 256 x 128 lattice points comprises 2nc particles. Due to lattice 
symmetry, it suffices to solve the problem in a single elementary cell. For p{r, t = 0) we 



choose a superposition of sharply localized Gaussians according to equation ( 13 ) with a 
large localization strength of a* = 200. 

At first, we fix n^ = 1 and calculate the equilibrium density profiles and the 
according approximate Helmholtz free energies for various coupling constants < 
F < 62.5. In figure 10 we plot the difference in Helmholtz free energy density 
AF(F)/A^ = F[p{r,t -^ oo;F)]/A^ - /o(F) between the final (solid/hquid) and the 
liquid state as a function of F. The system remains crystalline for couplings F > 30.7. 
However, the free energy difference is negative only for F > 36.2, which is equivalent 
with thermodynamic stability. As for the Gaussian parametrization coexistence is found 
in a narrow gap around F ^ 36.2 which we do not specify here. 

In figure [llta) we plot the equilibrium density profile p{r; F) for F = 36 which is 
close to freezing. In figure OHb) the quantity r'^p{r), where r is the distance from a 
lattice vector, is shown along the two directions [11] and [10], corresponding to cuts 
through the density plane in figure [ll|(a) along the x- and the y- axis, respectively, 
which is compared to a Gaussian of the same height as the density peaks. It is found 
that the density profile has an isotropic Gaussian form for small distances from the 
origin r < 0.1/p^/^. For larger distances, however, i.e., where the density is of the order 
p{r) < p, the density profile significantly deviates from a Gaussian form. In particular, 
we observe the establishment of "bridges" of higher density between neighbouring lattice 
sites, whereas the density is significantly lower between next-nearest neighbours. This 
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Figure 11. (a) The density profile p(r) obtained from dynamical DFT using 
the SOT with the RY closure for F = 36 which is close to freezing, (b) The 
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from the center to the respective edge of the box in (a) . The two curves are compared 
to a Gaussian of the same amplitude at r = 0. The inset displays the bare density 
along the same lines and the bare Gaussian. 



5 

C 
It 



0.04 
0.03 
0.02 
0.01 

-0.01 
-0.02 
-0.03 
-0.04 
-0.05 
-0.06 




r=35 
r=36 
r=37 
r=40 



n - 



0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 



Figure 12. The difference in Helmholtz free energy per particle AF{nc;T)/N as a 
function of ric for different coupling constants T = 35, 37, 40. The arrows indicate the 
positions of the minima. 



counter-intuitive behavior was also found applying the MWDA to hard sphere crystals in 
three spatial dimensions [H] . However, computer simulations revealed that the behavior 
should be the opposite. Although we did not measure the density profiles of the two- 
dimensional dipolar system in computer simulations, we expect a similar behavior: The 
probability density should be enhanced along the [ll]-direction as compared to the 
[10]-direction. 

We also performed the minimization procedure for different vacancy concentrations. 



In figure 12 we show the free energy difference AF(nc; T) as a function of Uc for four 
different values of F. We find, that for crystals in equilibrium, i.e., for F > 36, the 
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equilibrium vacancy concentration is 1 — ric ~ 0.03. However, the overheated crystal 
which is metastable for 31 < F < 36 prefers a vacancy concentration of 1 — n^ ~ —0.03, 
implying interstitials instead of vacancies. We note that most of the point defects in 
the experimental realization of the dipolar system appear in pairs or in pairs of pairs as 
dislocations or pairs of dislocations, respectively 



6. Discussion and concluding remarks 

In conclusion, we have demonstrated that density functional theory is able to 
quantitatively predict the freezing transition of a two-dimensional colloidal system 
with long-range 1/r ^-interactions in good agreement with experimental and simulation 
data. In complete analogy to systems in 3D, the appearance of long-range interactions 
requires the explicit inclusion of three-particle correlation functions of the liquid in the 
construction of the weighted density [IHl [22]. Furthermore, the predicted transition 
temperatures are very sensitive towards slight changes of the two- and three-particle 
correlation functions of the underlying fluid. A highly accurate input of the same is 
therefore crucial. 

The obtained density functional can be used in future studies in order to approach 
more complicated situations such as crystals in confinement [18], under gravity |17], 
and crystal- fluid interfaces [IHl HE]- By extending the static functional to Brownian 
dynamics [12], HHl HH HHl HE], one may even address nonequilibrium situations. One 
possible problem to tackle is heterogeneous nucleation upon temperature quenches and 
subsequent crystal growth as outlined recently in Ref. [IB] . 
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